第十章: 潛在成長模型分析

2024 三月 03

#整體設定,含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")

資料與管理

#讀取檔案
dta <- read.csv('../Data/KIT_memory.csv', 
                header=T, stringsAsFactors = TRUE)
#看資料結構
#程式報表10.1
dim(dta)
[1] 2000    7
#檢視資料格式與前六筆
head(dta)
識別碼 母親教育程度 月_03 月_06 月_12 月_18 月_24
S1001 大學以上 20 25 48 51 53
S1002 大學以上 25 29 37 50 59
S1003 專科以下 17 20 40 50 50
S1004 大學以上 25 28 46 50 59
S1005 大學以上 20 25 43 52 64
S1006 專科以下 19 27 42 46 57
#寬資料轉成長資料
#也另外製造「月」,以及將月置中的「月置中」變項
dtaL <- dta |> 
  tidyr::pivot_longer(cols=3:7,
                      names_to ='月齡',  
                      values_to = '記憶分數') |>
  dplyr::mutate(月 = parse_number(月齡), 
                月置中 = scale(月, scale = F)[,1])
#程式報表10.2
head(dtaL)
識別碼 母親教育程度 月齡 記憶分數 月置中
S1001 大學以上 月_03 20 3 -9.6
S1001 大學以上 月_06 25 6 -6.6
S1001 大學以上 月_12 48 12 -0.6
S1001 大學以上 月_18 51 18 5.4
S1001 大學以上 月_24 53 24 11.4
S1002 大學以上 月_03 25 3 -9.6
#看五波的平均數與變異數,母親教育程度的影響在月齡較大時顯現
#程式報表10.3
dtaL |> 
  tidyr::unite("月齡母親教育程度", 月齡, 母親教育程度) |>
  dplyr::select("月齡母親教育程度", "記憶分數") |>
  gtsummary::tbl_summary(by=月齡母親教育程度,
            digits=list(記憶分數 ~ c(1, 2)),
            type = all_continuous() ~ "continuous2",
            statistic = list(all_continuous() ~ "{mean} ({sd})")) |>                          
  modify_header(label ~ "**變項**") 
#看看偏態與峰度
#程式報表10.4前
dtaL |> 
  dplyr::group_by(月齡, 母親教育程度) |>
  dplyr::reframe(平均 = mean(記憶分數),
                 標準差 = sd(記憶分數),
                 偏度 = moments::skewness(記憶分數),
                 峰度 = moments::kurtosis(記憶分數))
月齡 母親教育程度 平均 標準差 偏度 峰度
月_03 大學以上 20.43 2.707 0.7297 3.019
月_03 專科以下 20.39 2.752 0.9074 3.579
月_06 大學以上 25.41 3.014 -0.1002 2.597
月_06 專科以下 25.42 3.133 -0.1345 2.617
月_12 大學以上 41.91 4.901 -0.6666 3.220
月_12 專科以下 41.12 5.481 -0.5547 4.134
月_18 大學以上 48.66 4.132 -0.1591 4.933
月_18 專科以下 47.21 4.565 -1.0957 9.287
月_24 大學以上 55.86 5.797 0.1587 2.332
月_24 專科以下 52.71 5.772 -0.1325 5.485
#歷次分數間相關
#程式報表10.4後
dta |> 
  dplyr::select(where(is.numeric)) |> 
  cor() |> 
  round(3)
月_03 月_06 月_12 月_18 月_24
月_03 1.000 0.483 0.328 0.252 0.177
月_06 0.483 1.000 0.438 0.351 0.248
月_12 0.328 0.438 1.000 0.525 0.396
月_18 0.252 0.351 0.525 1.000 0.579
月_24 0.177 0.248 0.396 0.579 1.000

繪圖

#繪製不同母親教育程度五波圖
#圖10.1a
p1 <- ggplot(data = dtaL, 
       aes(x = 月, y = 記憶分數)) +
 ggbeeswarm::geom_quasirandom(aes(group=月), alpha=.1)+
 facet_wrap(vars(母親教育程度))+
 stat_summary(fun = mean, geom = 'line', linewidth=.8) +
 scale_x_continuous(breaks=c(3,6,12,18,24))+
 labs(x = '月齡', 
      y = '記憶能力分數') +
 theme(legend.position = 'top')
#看五波的機率密度函數圖
#圖10.1b
p2 <- ggplot(dtaL, 
       aes(x=記憶分數, y=月齡)) +
  ggridges::geom_density_ridges2()+
  scale_x_continuous(breaks=seq(10, 70, by=10))+
  facet_wrap(vars(母親教育程度))+
  labs(x = '記憶能力分數', 
       y = '月齡')
p1+p2
#不同母親教育程度下的實徵累積機率密度圖
#圖10.2
ggplot(dtaL,
       aes(x=記憶分數, group=母親教育程度)) +
  stat_ecdf()+
  facet_wrap(vars(月齡))+
  labs(x="記憶能力分數",
       y="實徵累積機率密度")
#看看個別資料,畫上迴歸線並配上區間
#圖10.3
ggplot(data = dtaL, 
       aes(y=記憶分數, x=月)) +
 geom_line(aes(group = 識別碼), 
           linewidth = .2,
           alpha=.5,
           col=8) +
 stat_summary(fun.data = 'mean_cl_boot', geom = "errorbar", width = 0.5) +
 stat_smooth(method="lm", 
             formula= y ~ poly(x, 2),
             se=FALSE, 
             col=1, 
             linewidth=.5)+
 facet_wrap(vars(母親教育程度))+
 scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
 scale_y_continuous(breaks=seq(8, 72, by=8))+
 labs(x = '月齡', 
      y = '記憶能力分數')
#前述圖形看起來不完全是線性,後續進一步探索
#將資料框轉成 tsibble,方便後續使用
dta_tb <- brolgar::as_tsibble(x = dtaL,
                             key = 識別碼,
                             index = 月,
                             regular = FALSE)
#brolgar套件的features功能,可以計算每個幼兒記憶分數的五數綜合(含中位數)
#利用gghighlight,在繪圖時,特別強調中位數特別高或低的幼兒資料
#圖10.4
dta_tb |>
  brolgar::features(記憶分數, feat_five_num) |>
  dplyr::left_join(dta_tb, by = '識別碼') |>
  ggplot(aes(x = 月,
             y = 記憶分數,
             group = 識別碼)) +
  geom_line(linewidth=.2) + 
  gghighlight::gghighlight(med > 50 | med < 28, 
                           label_key = 識別碼,
                           use_direct_label = FALSE)+
   scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
  labs(x = '月齡', 
      y = '記憶能力分數')
#brolgar套件中的facet_strata(),預設將資料分成12組,助於探索資料
#看看個別資料
#圖10.5
dta_tb |>
  ggplot(aes(x = 月,
             y = 記憶分數,
             group = 識別碼)) + 
  #geom_point()+
  geom_line(alpha=.2, linewidth=.2)+
  brolgar::facet_strata()+
  scale_x_continuous(breaks=c(3, 6, 12, 18, 24))
#隨機取23位來看一下資料
#前面觀察到資料並非線性,嘗試以二次式迴歸配適資料
#為對比二次式,也用樣條迴歸(spline regression)配適資料
#圖10.6
set.seed(18072023)
dta |> dplyr::slice_sample(n=23) |>
  tidyr::pivot_longer(cols=3:7,
                      names_to ='月齡',  
                      values_to = '記憶分數') |>
  dplyr::mutate(月 = parse_number(月齡)) |>
  ggplot()+
  aes(x=月, y=記憶分數)+
  geom_point(pch=1, size=rel(3))+
  stat_smooth(method='lm', formula= y ~ poly(x, 2), 
              col=8,
              alpha=.5,
              linewidth=.6,
              se=FALSE)+
  stat_smooth(method='lm', 
         formula = y ~ splines::bs(x, knots = c(3, 12, 24), degree = 1),
              col=2,
              linetype='dotted',
              linewidth=1,
              se=FALSE)+
  facet_wrap(vars(識別碼)) +
  scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
  scale_y_continuous(breaks=seq(16, 72, by=16))+
  labs(x = '月齡', 
       y = '記憶能力分數',
       title="實灰:二次多項式, 點紅:樣條線性迴歸")

二次多項式迴歸

#以nlme,做每個人的迴歸線
m1 <- nlme::lmList(記憶分數 ~ 月置中 + I(月置中^2) | 識別碼, data = dtaL, subset=母親教育程度=='專科以下') |> 
  coef() |> 
  as.data.frame() |>
  dplyr::mutate(母親教育程度='專科以下')

m2 <- nlme::lmList(記憶分數 ~ 月置中 + I(月置中^2) | 識別碼, data = dtaL, subset=母親教育程度=='大學以上') |> 
  coef() |> 
  as.data.frame() |>
  dplyr::mutate(母親教育程度='大學以上')

m12 <- rbind(m1, m2)
names(m12)[1:3] <- c("b0", "b1", "b2")
#截距、一次項與二次項係數間的相關
cor(m12[,-4])
b0 b1 b2
b0 1.0000 0.4255 -0.6611
b1 0.4255 1.0000 0.0185
b2 -0.6611 0.0185 1.0000
#看看母親教育程度對係數的影響
#程式報表10.5
m12 |> 
  gtsummary::tbl_summary(by=母親教育程度,
                         statistic = list(all_continuous() ~ "{mean} ({sd})"))
Characteristic 大學以上, N = 1,3071 專科以下, N = 6931
b0 41.3 (3.9) 40.6 (4.3)
b1 1.80 (0.26) 1.67 (0.26)
b2 -0.05 (0.04) -0.05 (0.04)
1 Mean (SD)
#不同母親教育程度所得之截距,線性斜率,二次項效果估計值
#圖10.7
g1<-ggplot(m12,
       aes(x=b0, y=b1, color=母親教育程度))+
  geom_point(pch=1)+
  stat_ellipse()+
  scale_color_grey()+
  labs(x='截距估計值', 
       y = '線性斜率估計值')+
  theme(legend.position='top')
g2<-ggplot(m12,
       aes(x=b1, y=b2, color=母親教育程度))+
  geom_point(pch=1)+
  stat_ellipse()+
  scale_color_grey()+
  labs(x='線性斜率估計值', 
       y = '二次項效果估計值')+
  theme(legend.position='top')
g3 <- ggplot(m12,
       aes(x=b0, y=b2, color=母親教育程度))+
  geom_point(pch=1)+
  stat_ellipse()+
  scale_color_grey()+
  labs(x='截距估計值', 
       y = '二次項效果估計值')+
  theme(legend.position='top')
g1 + g2 / g3

潛在成長模型

先分析母親教育程度為大學以上者

#潛在成長模型可以視為因素分析的子模型,我們以 lavaan 進行分析
#先選擇母親教育程度為大學以上者分析
dta_2 <- dta |> dplyr::filter(母親教育程度 == '大學以上')

線性模型

#先試試看線性模型
growth1 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
slope =~ 0*月_03 + 1*月_06 + 3*月_12 + 5*月_18 + 7*月_24
'
#配適模型
fit1 <- lavaan::growth(model = growth1, data = dta_2)
#程式報表10.6前,後
fit1 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
   chisq       df   pvalue    rmsea     srmr      cfi      tli 
2039.065   10.000    0.000    0.394    0.352    0.000   -0.267 
fit1 |> parameterestimates() 
lhs op rhs est se z pvalue ci.lower ci.upper
intercept =~ 月_03 1.0000 0.0000 NA NA 1.0000 1.0000
intercept =~ 月_06 1.0000 0.0000 NA NA 1.0000 1.0000
intercept =~ 月_12 1.0000 0.0000 NA NA 1.0000 1.0000
intercept =~ 月_18 1.0000 0.0000 NA NA 1.0000 1.0000
intercept =~ 月_24 1.0000 0.0000 NA NA 1.0000 1.0000
slope =~ 月_03 0.0000 0.0000 NA NA 0.0000 0.0000
slope =~ 月_06 1.0000 0.0000 NA NA 1.0000 1.0000
slope =~ 月_12 3.0000 0.0000 NA NA 3.0000 3.0000
slope =~ 月_18 5.0000 0.0000 NA NA 5.0000 5.0000
slope =~ 月_24 7.0000 0.0000 NA NA 7.0000 7.0000
月_03 \~~ 月_03 3.1163 0.2533 12.301 0e+00 2.6198 3.6129
月_06 \~~ 月_06 5.5443 0.2739 20.241 0e+00 5.0075 6.0812
月_12 \~~ 月_12 40.5503 1.6438 24.669 0e+00 37.3285 43.7720
月_18 \~~ 月_18 7.0278 0.5132 13.694 0e+00 6.0219 8.0336
月_24 \~~ 月_24 31.4282 1.5183 20.700 0e+00 28.4524 34.4040
intercept \~~ intercept 4.1679 0.2834 14.705 0e+00 3.6124 4.7235
slope \~~ slope 0.3144 0.0261 12.024 0e+00 0.2631 0.3656
intercept \~~ slope -0.2282 0.0644 -3.544 4e-04 -0.3545 -0.1020
月_03 ~1 0.0000 0.0000 NA NA 0.0000 0.0000
月_06 ~1 0.0000 0.0000 NA NA 0.0000 0.0000
月_12 ~1 0.0000 0.0000 NA NA 0.0000 0.0000
月_18 ~1 0.0000 0.0000 NA NA 0.0000 0.0000
月_24 ~1 0.0000 0.0000 NA NA 0.0000 0.0000
intercept ~1 20.4577 0.0701 291.897 0e+00 20.3203 20.5950
slope ~1 5.5142 0.0214 257.357 0e+00 5.4722 5.5562

二次多項式模型

#試試看二次模型

growth2 <- '
intercept =~ 1*月_03 + 1*月_06    + 1*月_12 + 1*月_18 +  1*月_24
linear =~ (-3)*月_03 + (-2)*月_06 + 0*月_12 + 2*月_18 +  4*月_24
qudratic =~  9*月_03 + 4*月_06    + 0*月_12 + 4*月_18 + 16*月_24
'
fit2 <- growth(model = growth2, data = dta_2)
Warning in lav_object_post_check(object): lavaan WARNING: some estimated lv
variances are negative
#程式報表10.7前
fit2 |> parameterestimates() |>
  filter(!is.na(z)) 
lhs op rhs est se z pvalue ci.lower ci.upper
月_03 \~~ 月_03 3.4138 0.4015 8.5035 0.0000 2.6269 4.2006
月_06 \~~ 月_06 8.9889 0.4166 21.5760 0.0000 8.1723 9.8054
月_12 \~~ 月_12 23.1927 0.9902 23.4226 0.0000 21.2520 25.1335
月_18 \~~ 月_18 4.6782 0.4845 9.6562 0.0000 3.7286 5.6277
月_24 \~~ 月_24 22.6209 1.7112 13.2194 0.0000 19.2670 25.9748
intercept \~~ intercept 6.9909 0.5519 12.6672 0.0000 5.9092 8.0726
linear \~~ linear 0.3478 0.0260 13.4004 0.0000 0.2969 0.3986
qudratic \~~ qudratic -0.0139 0.0058 -2.3825 0.0172 -0.0253 -0.0025
intercept \~~ linear 1.1746 0.0998 11.7751 0.0000 0.9791 1.3702
intercept \~~ qudratic -0.0202 0.0424 -0.4769 0.6334 -0.1033 0.0629
linear \~~ qudratic -0.0420 0.0099 -4.2283 0.0000 -0.0615 -0.0226
intercept ~1 39.1876 0.0993 394.5548 0.0000 38.9929 39.3823
linear ~1 5.4303 0.0212 256.3939 0.0000 5.3888 5.4718
qudratic ~1 -0.3119 0.0084 -37.1707 0.0000 -0.3284 -0.2955
#修改二次模型,要求二次參數變異數為小的正數
#程式報表10.10後
growth21 <- '
intercept =~ 1*月_03 + 1*月_06    + 1*月_12 + 1*月_18 +  1*月_24
linear =~ (-3)*月_03 + (-2)*月_06 + 0*月_12 + 2*月_18 +  4*月_24
qudratic =~  9*月_03 + 4*月_06    + 0*月_12 + 4*月_18 + 16*月_24
qudratic~~.01*qudratic
'
fit21 <- growth(model = growth21, data = dta_2)
#程式報表10.7中,後
fit21 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
   chisq       df   pvalue    rmsea     srmr      cfi      tli 
1191.427    7.000    0.000    0.360    0.217    0.260   -0.057 
fit21 |> parameterestimates() |>
  filter(!is.na(z)) 
lhs op rhs est se z pvalue ci.lower ci.upper
月_03 \~~ 月_03 3.0453 0.3832 7.947 0.0000 2.2942 3.7963
月_06 \~~ 月_06 8.8895 0.4124 21.555 0.0000 8.0812 9.6977
月_12 \~~ 月_12 22.4378 0.9550 23.496 0.0000 20.5661 24.3095
月_18 \~~ 月_18 4.8065 0.4760 10.099 0.0000 3.8737 5.7394
月_24 \~~ 月_24 17.6857 1.2196 14.501 0.0000 15.2953 20.0761
intercept \~~ intercept 7.7949 0.5302 14.703 0.0000 6.7559 8.8340
linear \~~ linear 0.3694 0.0250 14.761 0.0000 0.3204 0.4185
intercept \~~ linear 1.0969 0.0986 11.127 0.0000 0.9037 1.2901
intercept \~~ qudratic -0.1533 0.0287 -5.339 0.0000 -0.2096 -0.0970
linear \~~ qudratic -0.0300 0.0100 -2.998 0.0027 -0.0496 -0.0104
intercept ~1 39.2071 0.1015 386.430 0.0000 39.0082 39.4060
linear ~1 5.4253 0.0211 256.864 0.0000 5.3839 5.4667
qudratic ~1 -0.3126 0.0088 -35.390 0.0000 -0.3299 -0.2953
#線性模型與二次模型比較
#程式報表10.8
anova(fit1, fit21)
Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(\>Chisq)
fit21 7 36091 36158 1191 NA NA NA NA
fit1 10 36933 36984 2039 847.6 0.4641 3 0

包含形狀因素不固定的潛在成長模型

#包含形狀因素的潛在成長模型,形狀不固定
growth3 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
shape =~     0*月_03 + 月_06   +   月_12 + 月_18   + 1*月_24
'
fit3 <- growth(model = growth3, data = dta_2)
#程式報表10.9前,後
fit3 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
 chisq     df pvalue  rmsea   srmr    cfi    tli 
60.079  7.000  0.000  0.076  0.045  0.967  0.953 
fit3 |> parameterestimates() |>
  filter(!is.na(z)) 
lhs op rhs est se z pvalue ci.lower ci.upper
shape =~ 月_06 0.1405 0.0022 65.114 0e+00 0.1363 0.1448
shape =~ 月_12 0.6066 0.0036 168.957 0e+00 0.5995 0.6136
shape =~ 月_18 0.7974 0.0032 245.741 0e+00 0.7910 0.8038
月_03 \~~ 月_03 3.1535 0.2450 12.873 0e+00 2.6734 3.6337
月_06 \~~ 月_06 4.9861 0.2506 19.893 0e+00 4.4948 5.4773
月_12 \~~ 月_12 13.9670 0.6157 22.686 0e+00 12.7604 15.1737
月_18 \~~ 月_18 5.4090 0.3997 13.532 0e+00 4.6256 6.1925
月_24 \~~ 月_24 18.5165 0.9268 19.979 0e+00 16.7000 20.3330
intercept \~~ intercept 4.2694 0.2780 15.356 0e+00 3.7245 4.8143
shape \~~ shape 15.2976 1.0050 15.222 0e+00 13.3279 17.2674
intercept \~~ shape -1.4559 0.3996 -3.644 3e-04 -2.2391 -0.6728
intercept ~1 20.4325 0.0753 271.236 0e+00 20.2848 20.5801
shape ~1 35.4099 0.1681 210.690 0e+00 35.0805 35.7393
#程式報表10.10
semTools::compareFit(linear = fit1, quadratic = fit21, shape = fit3, nested=F) |> summary()
####################### Model Fit Indices ###########################
              chisq df pvalue rmsea   cfi     tli  srmr        aic        bic
quadratic 1191.427   7   .000 .360  .260  -0.057  .217  36090.998  36158.280 
shape       60.079†  7   .000 .076† .967†   .953† .045† 34959.650† 35026.931†
linear    2039.065  10   .000 .394  .000  -0.267  .352  36932.636  36984.391 
#繪製路徑圖
#將圖存為 .png
#圖10.8
pl_01 <- lavaanPlot(model=fit3,
           edge_options = list(color = "grey"),
           coefs = TRUE,
           stand = FALSE,
           covs = TRUE)
dir.create(file.path(here::here()), 'figure')
lavaanPlot::save_png(pl_01, "figure/plot_01.png")
knitr::include_graphics("figure/plot_01.png")

包括共變量的潛在成長模型

#我們選擇包含形狀因素的潛在成長模型
#進一步看看母親教育程度對截距與形狀的影響。
#先看看不同母親教育程度下,資料配適情形。
dta_1 <- dta |> 
  dplyr::filter(母親教育程度 == '專科以下')
fit32 <- growth(model = growth3, data=dta_1)
#程式報表10.11前,後
fit32 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
 chisq     df pvalue  rmsea   srmr    cfi    tli 
24.619  7.000  0.001  0.060  0.041  0.980  0.972 
fit32 |> parameterestimates() |>
  filter(!is.na(z)) 
lhs op rhs est se z pvalue ci.lower ci.upper
shape =~ 月_06 0.1559 0.0033 47.675 0.0000 0.1495 0.1624
shape =~ 月_12 0.6417 0.0059 109.625 0.0000 0.6302 0.6531
shape =~ 月_18 0.8302 0.0050 165.500 0.0000 0.8203 0.8400
月_03 \~~ 月_03 3.2822 0.3500 9.379 0.0000 2.5963 3.9681
月_06 \~~ 月_06 5.0796 0.3524 14.415 0.0000 4.3889 5.7703
月_12 \~~ 月_12 18.0393 1.0802 16.700 0.0000 15.9222 20.1565
月_18 \~~ 月_18 7.2380 0.6609 10.952 0.0000 5.9427 8.5334
月_24 \~~ 月_24 15.8757 1.1730 13.534 0.0000 13.5766 18.1748
intercept \~~ intercept 4.4440 0.3991 11.134 0.0000 3.6617 5.2262
shape \~~ shape 15.6778 1.4137 11.090 0.0000 12.9070 18.4485
intercept \~~ shape -0.8501 0.5604 -1.517 0.1293 -1.9485 0.2483
intercept ~1 20.3867 0.1055 193.177 0.0000 20.1799 20.5936
shape ~1 32.3128 0.2241 144.212 0.0000 31.8736 32.7520
#將母親教育程度變成 0,1 的虛擬變項
dta <- dta |> 
  dplyr::mutate(母親大學以上 = ifelse(母親教育程度 =='大學以上', 1 , 0))

包含形狀因素的潛在成長模型

#包含形狀因素的潛在成長模型,讓性別影響截距與形狀因素。
growth4 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
shape =~     0*月_03 + 月_06   + 月_12  + 月_18   + 1*月_24
intercept ~ 母親大學以上
shape ~ 母親大學以上
'
fit4 <- growth(model = growth4, data = dta, fixed.x = T)
#程式報表10.12前,後
fit4 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
  chisq      df  pvalue   rmsea    srmr     cfi     tli 
125.850  10.000   0.000   0.076   0.040   0.957   0.935 
fit4 |> parameterestimates() |>
  filter(!is.na(z)) 
lhs op rhs est se z pvalue ci.lower ci.upper
shape =~ 月_06 0.1453 0.0018 80.435 0.0000 0.1418 0.1489
shape =~ 月_12 0.6176 0.0031 199.605 0.0000 0.6115 0.6236
shape =~ 月_18 0.8077 0.0027 293.860 0.0000 0.8023 0.8130
intercept ~ 母親大學以上 -0.1331 0.1202 -1.107 0.2683 -0.3687 0.1026
shape ~ 母親大學以上 2.2623 0.2330 9.709 0.0000 1.8056 2.7190
月_03 \~~ 月_03 3.2297 0.2015 16.027 0.0000 2.8348 3.6247
月_06 \~~ 月_06 5.0261 0.2047 24.552 0.0000 4.6249 5.4273
月_12 \~~ 月_12 15.4539 0.5485 28.174 0.0000 14.3789 16.5290
月_18 \~~ 月_18 6.0735 0.3485 17.430 0.0000 5.3906 6.7565
月_24 \~~ 月_24 17.9492 0.7409 24.227 0.0000 16.4971 19.4013
intercept \~~ intercept 4.3130 0.2282 18.897 0.0000 3.8657 4.7604
shape \~~ shape 15.3526 0.8222 18.672 0.0000 13.7411 16.9642
intercept \~~ shape -1.2111 0.3257 -3.718 0.0002 -1.8494 -0.5727
intercept ~1 20.5061 0.0995 206.120 0.0000 20.3111 20.7011
shape ~1 32.8697 0.2023 162.507 0.0000 32.4733 33.2662

比較實際記憶分數平均與預測值平均

#比較實際資料與預測值平均
#計算實際平均數
df <- aggregate(cbind(月_03, 月_06, 月_12, 月_18, 月_24) ~ 母親教育程度, 
                data = dta, mean) |>
  tidyr::pivot_longer(cols=-1, 
                      names_to = "月", 
                      values_to = '記憶分數平均') |>
  dplyr::mutate(月 = parse_number(月))
#計算預測值
t <- c(0,.145, .618, .808, 1)
df$估計值平均 <- c((20.506-.133)+(32.870+2.262)*t, 20.506+32.870*t)
#畫圖
#圖10.9
df |> tidyr::pivot_longer(cols=3:4,
                      names_to = "估計",
                      values_to = "平均分數") |>
 ggplot()+
   aes(y=平均分數, x=月, group=估計, color=估計)+
   geom_point(size=rel(5), pch=1)+
   geom_line(linewidth=.4)+
   facet_wrap(vars(母親教育程度))+
   scale_color_grey(end=.6)+
   scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
   labs(x = '月齡',
        y = '記憶能力平均分數')+
   theme(legend.position='top')
#繪製路徑圖
#將圖存為 .png
#圖10.10
pl_02 <- lavaanPlot(model=fit4,
           edge_options = list(color = "grey"),
           coefs = TRUE,
           stand = FALSE,
           covs = TRUE)
dir.create(file.path(here::here()), 'figure')
lavaanPlot::save_png(pl_02, "figure/plot_02.png")
knitr::include_graphics("figure/plot_02.png")